
/*
 * GA.d: Genetic Algorithms and helpers.
 * Author: Guilherme R. Lampert
 * 2013-11-25
 */

import std.stdio;
import std.random;
import std.bitmanip;
import std.algorithm;

// Standard values used by the GA:
enum NUM_ELITE = 4;
enum NUM_COPIES_ELITE = 1;

int BinToInt(in ubyte[] bits)
{
	int val = 0;
	int multiplier = 1;

	for (int b = cast(int)bits.length; b > 0; --b)
	{
		val += (bits[b - 1] * multiplier);
		multiplier *= 2;
	}

	return (val);
}

double RandFloat()
{
	// Random float between 0 and 1:
	return uniform!("[]", double, double)(0.0, 1.0);
}

// ======================================
// Genome:
// ======================================

class Genome {

	// Genome bits (each byte is treated as a bit (0|1)):
	ubyte[] bits;

	// Fitness value, from 0 to 1.
	double fitness;

	// How many bits per gene:
	int geneLength;

	// Construct with random gene bits:
	this(int numBits, int geneSize)
	{
		this.fitness = 0.0;
		this.geneLength = geneSize;

		// Fill with random bits:
		this.bits = new ubyte[numBits];
		for (int i = 0; i < numBits; ++i)
		{
			this.bits[i] = uniform!("[]", ubyte, ubyte)(0, 1); // random between 0 and 1
		}
	}

	// Construct with provided gene bits:
	this(int numBits, int geneSize, ubyte[] chromosomeBits)
	{
		assert(numBits == chromosomeBits.length);

		this.fitness = 0.0;
		this.geneLength = geneSize;
		this.bits = chromosomeBits;
	}

	// Decodes each gene into decimal:
	int[] Decode()
	{
		int[] decoded = new int[bits.length / geneLength];
		int i = cast(int)(decoded.length - 1);

		// Step through the chromosome a gene at a time:
		for (int gene = 0; gene < bits.length; gene += geneLength)
		{
			// Get the gene at this position:
			ubyte[] thisGene = bits[gene .. (gene + geneLength)];

			// Convert to decimal and add to list of decoded:
			assert(i >= 0);
			decoded[i--] = BinToInt(thisGene);
		}

		return (decoded);
	}

	// Encode data into the gene's binary string:
	void Encode(in int[] data)
	{
		ubyte * x = cast(ubyte *)data.ptr;
		ulong nBits = (data.length * int.sizeof * 8); // total number of bits.
		ulong shift = 0;

		assert(data.length == (bits.length / geneLength));
		assert(nBits == bits.length);

		while (nBits--)
		{
			bits[nBits] = (((*x >> shift++) & 1) ? 1 : 0);
			if (shift >= 8)
			{
				++x;
				shift = 0;
			}
		}
	}

	// Comparison operator for sorting:
	int opCmp(ref const Genome genome) const
	{
		if (this.fitness > genome.fitness)
		{
			return (+1);
		}
		if (this.fitness < genome.fitness)
		{
			return (-1);
		}
		return (0); // Equal
	}
}

// ======================================
// GeneticAlgorithm:
// ======================================

class GeneticAlgorithm {

	// Data:
	Genome[] genomes;           // Population of genomes
	int      populationSize;    // Current population
	double   crossoverRate;     // Crossover rate
	double   mutationRate;      // Mutation rate
	int      chromoLength;      // How many bits per chromosome
	int      geneLength;        // How many bits per gene
	double   totalFitnessScore; // Fitness of all individuals
	int      generationNum;     // Generation number

	// Default constructor, with params:
	this(double crossRate, double mutRate, int popSize, int numBits, int geneLen)
	{
		this.populationSize    = popSize;
		this.crossoverRate     = crossRate;
		this.mutationRate      = mutRate;
		this.chromoLength      = numBits;
		this.geneLength        = geneLen;
		this.totalFitnessScore = 0;
		this.generationNum     = 0;

		CreateInitialPopulation(); // Create random population
	}

	// Default constructor, with initial population:
	this(double crossRate, double mutRate, int popSize, int numBits, int geneLen, Genome[] population)
	{
		this.populationSize    = popSize;
		this.crossoverRate     = crossRate;
		this.mutationRate      = mutRate;
		this.chromoLength      = numBits;
		this.geneLength        = geneLen;
		this.totalFitnessScore = 0;
		this.generationNum     = 0;
		this.genomes           = population;
	}

	// Iterates through each genome flipping the bits according to the mutation rate.
	void Mutate(ref ubyte[] bits)
	{
		for (int b = 0; b < bits.length; b++)
		{
			// Do we flip this bit?
			if (RandFloat() < mutationRate)
			{
				// Flip it:
				bits[b] = !bits[b];
			}
		}
	}

	// Takes 2 parent bit vectors, selects a midpoint and then swaps the ends
	// of each genome creating 2 new genomes which are stored in baby1 and baby2.
	void Crossover(in ubyte[] mum, in ubyte[] dad, out ubyte[] baby1, out ubyte[] baby2)
	{
		// Just return parents as offspring dependent on the rate or if parents are the same:
		if ((RandFloat() > crossoverRate) || (mum == dad))
		{
			baby1 = mum.dup;
			baby2 = dad.dup;
			return;
		}

		// Resize:
		baby1 = new ubyte[mum.length];
		baby2 = new ubyte[dad.length];

		// Determine a crossover point:
		const int cp = uniform(0, chromoLength);

		// Swap the bits:
		baby1[0 .. cp] = mum[0 .. cp];
		baby2[0 .. cp] = dad[0 .. cp];
		baby1[cp .. baby1.length] = dad[cp .. dad.length];
		baby2[cp .. baby2.length] = mum[cp .. mum.length];
	}

	// Selects a member of the population by using roulette wheel selection.
	Genome RouletteWheelSelection()
	{
		const double slice = (RandFloat() * totalFitnessScore);
		int	selectedGenome = 0;
		double total = 0.0;

		for (int i = 0; i < populationSize; ++i)
		{
			total += genomes[i].fitness;

			if (total > slice)
			{
				selectedGenome = i;
				break;
			}
		}

		return (genomes[selectedGenome]);
	}

	// This works like an advanced form of elitism by inserting 'numCopies'
	// copies of the NBest most fittest genomes into a population vector.
	void GrabNBest(int NBest, int numCopies, out Genome[] newPopulation)
	{
		sort(genomes);

		// Add the required amount of copies of the N fittest to the supplied vector:
		while (NBest--)
		{
			for (int i = 0; i < numCopies; ++i)
			{
				newPopulation ~= genomes[NBest];
			}
		}
	}

	// Creates the initial population:
	void CreateInitialPopulation()
	{
		// Clear existing population, if any
		genomes = new Genome[populationSize];

		for (int i = 0; i < populationSize; ++i)
		{
			genomes[i] = new Genome(chromoLength, geneLength);
		}

		// Reset all variables:
		generationNum = 0;
		totalFitnessScore = 0;
	}

	// Total fitness for population:
	void CalculateTotalFitness()
	{
		totalFitnessScore = 0;
		for (int i = 0; i < populationSize; ++i)
		{
			totalFitnessScore += genomes[i].fitness;
		}
	}

	// This is the workhorse of the GA. It first updates the fitness
	// scores of the population then creates a new population of
	// genomes using the Selection, Crossover and Mutation operators.
	bool Run()
	{
		Genome[] babyGenomes = [];
		CalculateTotalFitness();

		// Now to add a little elitism we shall add in some copies of the fittest genomes:
		// Make sure we add an EVEN number or the roulette wheel sampling will crash!
		if (!(NUM_COPIES_ELITE * NUM_ELITE % 2))
		{
			GrabNBest(NUM_ELITE, NUM_COPIES_ELITE, babyGenomes);
		}
		else
		{
			return (false);
		}

		for (;;)
		{
			// Select 2 parents:
			Genome mum = RouletteWheelSelection();
			Genome dad = RouletteWheelSelection();

			// Crossover:
			Genome baby1 = new Genome(chromoLength, geneLength);
			Genome baby2 = new Genome(chromoLength, geneLength);
			Crossover(mum.bits, dad.bits, baby1.bits, baby2.bits);

			// Mutate:
			Mutate(baby1.bits);
			Mutate(baby2.bits);

			// Append to new population:
			babyGenomes ~= baby1;
			if (babyGenomes.length >= populationSize) { break; }
			babyGenomes ~= baby2;
			if (babyGenomes.length >= populationSize) { break; }
		}

		// Copy babies back into starter population:
		genomes = babyGenomes;

		// Increment the generation counter:
		++generationNum;

		return (true);
	}
}
